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We present a numerical approach which allows the solving of Bethe equations whose solutions 
define the eigenstates of Gaudin models. By focusing on a new set of variables, the canceling di- 
vergences which occur for certain values of the coupling strength no longer appear explicitly. The 
problem is thus reduced to a set of quadratic algebraic equations. The required inverse transfor- 
mation can then be realized using only linear operations and a standard polynomial root finding 
algorithm. The method is applied to Richardson's fermionic pairing model, the central spin model 
and generalized Dicke model. 



The correspondence between Bethe Ansatz / Ordinary 
Differential Equation (ODE) has been found some time 
ago^ on the basis of similarity between T-Q system of the 
Bethe Ansatz and functional relations including spectral 
determinants on the ODE side. For extensive review we 
refer to^. The correspondence is related to the Langlands 
correspondence where the ODE part has to do with so- 
called Miura opers^. It has many interesting links with 
conformal field theory, quasi-exact solvability, supersym- 
metry and PT-symmetric quantum mechanics. It was 
also used it in the context of physics of cold atoms and 
quantum impurity problem"*. Here we develop further 
this remarkable correspondence to facilitate the solving 
of otherwise difHcult-to-solve Gaudin systems and imple- 
ment it to several physically interesting sitiiations. 

Indeed, numerous integrable models derived from a 
generalized Gaudin algebra have been used to describe 
the properties of physical systems. The Richardson 
fermionic pairing Hamiltonian^ has explained most prop- 
erties of superconducting nanograins^''', numerous stud- 
ies of the dccoherence of a single electron spin trapped in 
a quantum dot rely on the central spin model^ while the 
inhomogeneous Dick; models^ have been used to study 
light-matter interaction in many cavity-based systems. 

The quantum integrability of such systems brings ma- 
jor simplifications to the structure of their eigenstates. 
In a system containing iV degrees of freedom, they can 
be fully described by a number M of complex parameters 
we call rapidities. M is here a number of the same order 
as A'' despite the exponentially large Hilbert space. The 
possible values of this set of parameters can be obtained 
by finding sohitions to a set of M coupled non-linear al- 
gebraic equation: the Bethe equations. However, analyt- 
ically finding solutions to the Bethe equations remains 
impossible except in very specific cases. The problem 
is still numerically approachable but solving non-linear 
equations remains a challenge that only iterative meth- 
ods can tackle. 

While previous efforts in designing algorithms have 
shown promising results [10-14], methods which are suf- 
ficiently fast and stable for systematically finding a large 
number of eigenstates remained elusive. As was shown 



by one of the authors, computing efficiently only a small 
fraction of the complete Hilbert space can be sufficient 
to access static^^, dynamical*^ and even non-equilibrium 
dynamical properties of these systems*'^'** making such 
a method highly desirable. 

This paper presents an algorithm which finds solutions 
to the system of equations with unprecedented speed and 
stability. The first section discusses general consider- 
ations related to the quantum systems discussed here. 
Section II describes in detail the method used to solve 
for a set of intermediate variables and we show in section 
III, how one can recover the rapidities from these vari- 
ables. Explicit numerical applications to the Richardson 
model, the generalized Dicke model and the central spin 
model are presented in section IV and we conclude in the 
remaining section. 



I. THE MODELS 

Let us first introduce the generalized Gaudin algebra 
defined by the operators S^(Ai), S^(Ai), S^(Aj) with A^ 
any complex number and the commutation rules they 
obeyi9.20. 



[S-(AO,SJ'(Aj)] = i{Y{Xi,Xj)S'{Xi) - X{Xi,Xj)S'{Xj)), 

[S^(Ai),S^(A,-)] = i(^(A„A,)S^(A,)-y(A„A,)S^(A,)), 

[S^(A,),S^(A,)] = *(A(A„A,)S«(Ai)-^(Ai,Aj)S^(A,)), 

[S'^(Ai),S''(A,)] =0 K = x,y,z. (1) 

In this paper, we deal with the rational family of 
Gaudin models for which 



X{Xi,Xj)=Y{X„Xj) = Z{Xi,Xj) 



Xi — Xj 



■ (2) 



The generalized Dicke model we also consider is how- 
ever obtained from a trigonometric Gaudin model^* de- 
fined by 
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where 



X{X,,Xj) - y(A,,A,)-— ^ 
Z{Xi,Xj) = gcot{Xi-Xj). (3) 

The model is then derived by taking the large spin 

limit of a Holstcin-Primakoff transformation for one of 
the degrees of freedom. The resulting limit shares more 
resemblance to the rational family of Gaudin models and 
simply constitutes its extension to include a bosonic de- 
gree of freedom. We will therefore indiscriminately use 
the term rational to also include models derived from this 
particular limit. 

For any given realization of the Gaudin algebra, one 
can define a set of N commuting operators Ri allowing 
one to build exactly-solvable Hamiltonians 



N 



(4) 



M 



P{z) = \[{z-Xk) 



(8) 



fe=i 



is the polynomial of degree M whose M roots correspond 
to the values of A^. 

Since h.{z) obeys the following Riccati-type differential 
equation^^'^^ 



dK{z) 
dz 



+ K\z) 



+E 



1 



^^-^o^y - ^oc){z - Xp) 
y ' 

(7- 



{z - Xa){Xa - Xp) ' 



(9) 



it is easy to show that when the set {Aj} is a solution of 
the Bethe equations we have 



for which the Ri are evidently constants of motion. 
Eigcnstates of these models are, for a given number of 
excitations M, all defined by the construction 



|{Ai...AM})(xnS+(A.)|0). 



(5) 



Here S"'"(u) = S^(u) -|- zS^(m) is a family of generalized 
creation operators parametrized by the complex parame- 
ter u. Its explicit expression in terms of the fundamental 
operators defining a particular realization will be model 
dependent. The pseudovacuum |0) is defined as the low- 
est weight vector, i.e. S~{u) |0) = 0, and can also differ 
in distinct realizations of the Gaudin algebra. 

States of the form (5) become eigenstates of the system 
provided the M rapidities Aj are solution of a set of cou- 
pled non-linear algebraic equation: the Bethe equations. 
For rational models, these equations can be written, in 
general, as 



Fi^i) = E 



1 



Xi — A, 



(6) 



with S''{Xi)\0) = i^(Ai)|0) defining the lowest weight 
function F. 

One of the major difficulties in solving these equations 
numerically is the divergences which occur whenever two 
rapidities coincide. While they are cancelled by similarly 
diverging terms on the left hand side of the equation, 
it still has an important impact on numerical stability 
and computational speed. To circumvent these potential 
pitfalls, we introduce the function 



A(^) 



M 

E 

k=l 



(7) 



A'iz)+A^z)-J2^^ = 0. (10) 

One can derive this last equation with respect to z any 
number of times to write additional equalities: 

A"(.)+2A(z)A'(.) + ^^^^^=0, 



A"'{z) + 2A{z)A"{z) + 2A'{zf - 



(^-A„)3 



= 0, 



(11) 



From this point on, we will assume a particular form 
for the F function which, while not general, encompasses 
a wide variety^" of physically relevant realizations of the 
Gaudin algebra. We restrict ourselves to the following 



N 



A. 



+ —X+ — 
(ei-Aa) 2g " 2g' 



(12) 



The exact physical nature of the parameters g and 
Ci is highly model dependent, but in the cases treated 
in this paper it will be made explicit in section IV 
when discussing their specifics. In (pseudo-)spin models, 
Ai = \si\fli with I Si I the norm of the local spin degree of 
freedom, while is an integer relatable to the degener- 
acy, i.e. the number of elements of the set {ej} equal to 

. Consequently, every Ai can then take any integer or 
half-integer value. 



II. SOLVING THE SYSTEM 

From the previously found set of differential equations, 
we can write a new set of algebraic equations by simply 
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taking the limits z — > of Eqs. (13) and (11). Using the 
previous form of -F(Aa) (Eq. (12)) we find 



(l-2A,)A'(e,) + A2(e,) + -M 



A(e,-) 



A(e,)-A(e,) 



Ei - e,- 



(13) 



B 



(l-A,)A"(6,) + 2A(e,)A'(e,)--A(e,) 
_Be^ ^ A(e^2^A(eO 



2A,: 



(14) 



1 - A'"(e,) + 2A(e,)A"(e,-) + 2A'(e,)2 

-2|A'(e,) - ^^A"(.,) 



2A,A"(£,) 
e,; - 



^ "^J 



= 0, 



(15) 



One can immediately see that for any non-degenerate 

spin 1/2 degree of freedom {Aj = the first differen- 
tial equation reduces to a quadratic algebraic one which 
depends only on the set of variables {A{ej)}. This is due 
to the canceling of the first term (1 — 2Aj = 0). 

For Aj = 1, which can occur either due to a doubly 
degenerate spin 1/2 or a single spin 1, the two first equa- 
tions form a quadratic system of equations depending on 
{A{ej)} but also on the first derivative A'(ej) evaluated 
at ej. Larger values of Aj require additional equations, 
but it is always possible to write a closed coupled sys- 
tem of quadratic equations. It would depend on A(ej) 

on both A(ej), A'(ej) for vari- 



for variables with Aj = ^ 



ablcs with Aj = 1, A(ej), A'(ej), A"(ej) for variables with 
Aj ■ = I , etc. The resulting system of equations is built by 
using the n first equations above for any variables with 
Aj = I and its solutions would give a one to one corre- 
spondence with the solutions of the Bethe equations and 
therefore the eigenstates of the system. 

There is one caveat that readers should be aware of. 
While a spin S* > ^ or a set of 25' degenerate spin-^ lead 
to the same set of Bethe equations, only in the former 
case would the solutions to the Bethe equations give us 
the full Hilbert space. This is particularly simple to un- 
derstand for a system containing only one spin S = 1 



or two degenerate spin-^. In the first case, the Hilbert 
space dimension is 3 while in the second it is 4 and the 
resulting Bethe equations have only 3 distinct solutions. 
In the degenerate case, the Bethe equations would only 
give us highest weight states (one can think of it as the 
J = 1 triplet for two spins- 5) and the remainder of the 
Hilbert space would need to be reconstructed by building 
the appropriate set of orthogonal states. 

While in principle the system of quadratic equations 
can be solved for any spin or degeneracies, for the re- 
mainder of this paper will focus on the simplest case of 
non-degenerate spin 1/2 systems (Aj = | V j). In this 
case, the closed set of algebraic equations is given by the 
A'' following equations: 



f A(e,)-A(e.)_^B 

f^j ^j-^i 9 \9 ' 9 J 



A(e,). 
(16) 



In the subspace M < N this system of equations is 

larger than the original one (Eq. (6)). However, being 
quadratic and not having canceling divergences, makes it 
a much easier problem to tackle numerically. 

We also note in passing that the values of A(ej) de- 
termine the eigenvalues Vj of the commuting Hamil- 
tonians Rj. The eigenvalues of the transfer matrix, 
r(A) = |Tr(S^) are given in terms of F{\) (see Eq.(12)) 
as follows: 



r(A) = F\\) + F'(A) + 2 ^(^^(A) - F(A,))/(A - A,), 

(17) 

to {BM/2g) + 
eigenvalue Vj of 
the conserved operator Rj is given by the pole of t(A) 
at A = ej we can read off this eigenvalue by looking at 
the residue of r(A) at this pole 

r, = 2AjA{ej) - A,(C + B)/g + 2^ A,^,/(e, - e,). 



with the last term reducing 
X;,- AjA(ej)/(A - Cj). Since the 



(18) 



As with any non- linear system, solving Eqs. (16) ne- 
cessitates an iterative method such as the well known 
Newton-Raphson approach. Provided we have access to 
an initial approximative solution that is good enough to 
be in its basin of attraction, convergence of the method 
to a specific solution will be quadratic . Finding eigen- 
states therefore requires a sufficiently good guess which 
we obtain by slightly deforming known solutions. In the 
cases treated here, we exclusively know the exact eigen- 
states of the system at g = 0. We therefore approach 
the problem by slowly deforming these <? = solutions. 
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Provided g is raised in small enough steps, this guaran- 
tees that the previously found solutions can be used to 
generate a good approximation at the next point. 

Remarkably, it is numerically simple to compute the n 
first derivatives of the variables Aj — gA{ej) with respect 
to g. Indeed, one can show that every order needs only 
the solving of the same linear system with an updated 
right hand side that is straightforwardly computed from 
the previously computed information. Defining A^"^ = 
we obtain in matrix form the following linear system 



(19) 



for some right hand side that depends only on the lower 
derivatives of A(ej) and a constant matrix 

77^ 



The components of the right hand side Rn,j can be com- 
puted iteratively as 



Ri 



'0,j 



= -A? 



Rnj = ^(-i?„-i,,+Af-i)[(Be,+C)-2A,]) 



which gives expressions for the expectation values and 
correlation functions in terms of simple determinants^^ 
built out using the set of rapidities {A,;}. It is therefore 
necessary to be able to extract them from a given solution 
obtained in terms of the variables {A(ei)}. 

Going back to the definition of A(z) (Eq. 7), New- 
ton's identities allow us to write it explicitly in terms of 
elementary symmetric polynomials. We have 



M 



A(e, 



J2 mef-^Pi 

=0 

M 



M-m 



m=0 
M 



(23) 



M-m 



m=0 



with 



Pk 



E 



'^<h<h<-<3k<M 



Pn = 1. 



(24) 



Having the TV values A(ej) at hand, writing a lin- 
ear problem for the elementary symmetric polynomials 
is then a trivial matter: 



n 



(21) 

Since this requires a single matrix inversion (LU or QM 

decomposition to be exact) which can be reused at every 
order, it is a numerically a fast process to compute the n 
first derivatives. Using the resulting derivatives, a Tay- 
lor expansion gives an excellent initial approximation to 
A(ej) at g + Ag even for fairly large Ag: 



A(e,) 



A(e,)L 



^ A;! 

fe=i 



1 d^Aie,) 



{dgY 



{Agf. (22) 



In principle, the radius of convergence of the Taylor 
expansion around the current solution would set an up- 
per limit on the Ag step we can take. Nonetheless, one 
should keep in mind that adding terms to the Taylor se- 
ries only offers linear convergence while Newton's method 
converges quadratically. While computing less deriva- 
tives and relying on more Newton steps could speed up 
the calculation, it also comes with the risk that the ini- 
tial guess falls outside the basin of attraction of the de- 
sired solution. Optimizing the computation speed with- 
out compromising stability is then a question of balance 
between the number of computed derivatives and the size 
of the steps one takes in g. 



III. INVERTING A(z) 

Our capacity to compute physically relevant quanti- 
ties relies on the use of Slavnov's determinant formula^"* 



M-l 



^ {meJ-'-A{ei)ef)PM-m = A(e,)ef - Mej 



m=0 



(25) 



Case M < N 



When the number of excitations M is smaller than A'', 
one simply needs to pick M values of (cj) to extract the 
polynomials PM-m- The complete set of M elementary 
symmetric polynomials give the real coefficients of the 
single variable polynomial P{z) at every order. The cor- 
responding ensemble of parameters {A^} which defines 
the eigenstate is then obtained by finding every roots of 
P{z). 

This is a well studied problem for which many meth- 
ods exists. Here we choose to use Laguerre's method 
with polynomial deflation. Laguerre's method is an al- 
most "sure-fire" method for finding one root ri. We then 
proceed by finding a root r2 of the deflated polynomial 



at the next step, then a root of 7- 



etc. 



-ri)(z-r2)' 

Repeating M times allows the extraction of every root of 
our initial polynomial, i.e. every rapidity A^. 

This procedure for inverting the A(z) function is purely 
local in g in the sense that it only needs A(ej) at a given 
g. Therefore it is not necessary to perform the inversion 
at values of g we are not interested in. Solving for A{ej) 
still requires this scan in g but if we are, for example, only 
interested in the large g regime, not having to perform 
the inversion at intermediate steps can markedly reduce 
computation time. 
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B. Case M>N 

For a Hamiltonian realized in terms of only spin 1/2 
degrees of freedom the number of excitations is always 
M < N. However, in more general cases there exists a 
subspace where the excitation number is larger than the 
number of degrees of freedom (M > N). In such a case, 
one cannot simply use the set of N parameters A(ej) to 
find the M > N elementary symmetric polynomials. 

However, for any solutions with M < X^ili one 
can still extract them by solving only linear problems, 
since the supplementary information has been obtained 
from the necessity to solve for certain derivatives A'(ei). 
The unbounded number of excitations in models contain- 
ing bosonic degrees of freedom such as the Dickc mod- 
els presented here, would give a subset of states with 
M > X^^i 2^j. These would not allow such a simple 
inversion. However, since any degree of freedom with 
Ai can only accommodate up to 2Ai excitations, for any 
model based exclusively on spins (of any magnitude and 
degeneracy) the total number of excitations is always 
bounded and every solution has M < J2^=i 

In solving the quadratic equations (13,14,...) we are 
provided with A(ei) as well as its 2Ai — 1 first derivatives. 
Having these values at hand, writing a linear problem 
for the elementary symmetric polynomials is simply a 
matter of defining rational functions which are linear in 
Pn in both the numerator and denominator. Naturally, 
we already have from the definition 

N 

= ^ = ' (26) 

Z^PN-m 

m=0 

but we can also write 
V^{z)^K"{z) + 2,K{z)K'{z)+K{zf = 

(27) 

By evaluating the M necessary functions 
{A.[z),V\{z),V2{z), ..^) at z = Cj, we can write M 
equations linear in the polynomials P„. Solving them 
and using a root finding algorithm for P{z) would, once 
again, give us the set of rapidities {Aj}. 

IV. APPLICATIONS 

We now turn to three specific models derived from a 
generalized Gaudin algebra in order to demonstrate the 



efficiency of this approach. First we treat the discrete 
reduced BCS^*' model (Richardson model). One of the 
authors has already been involved in solving the Bethe 
equations in this context^^"^*. However, in this series of 
papers, it was performed without the currently discussed 
method. This led to serious stability and computation 
time issues which are now extremely well controlled in 
the current approach. 

A. Richardson model 

The Richardson model^, while having a rich history 
in nuclear physics, has also recently found renewed in- 
terest in condensed matter physics in light of tunneling 
experiments performed on superconducting nanograins^. 
It is nothing but a discrete version of the reduced ECS 
model which limits the interaction to a uniform s-wave 
pairing term between time reversed states. Using An- 
derson pseudospin-i representation in terms of fermionic 
operators 

0+ _ J J 
•^j — ''it 

Si = c^Cit, (28) 
the Hamiltonian is given by 

N N 

H = Y,^iSt-gY,StSj. (29) 

Here Cj corresponds to the discrete set of unblocked sin- 
gle fermion energies which can accommodate a Cooper 
pair, while g is the pairing strength between time- 
reversed fermionic states. While its integrability in a 
strict sense was proven much later^^, the exact solution 
of the model was first proposed by Richardson himself^. 
The eigenstates of the system are of the form given in 
Eq. (5) with a pseudovacuum |0) = ||,|, ...,|) that 
has fully down polarized pseudospins (Fock vacuum of 
Cooper pairs ) and the Gaudin creation operators given 
by 

N ^+ 

The Bethe equations for this particular reahzation are 
given by^^'^^ 

£a^ + ^^£^ . = 1,...,M. (31) 

This obviously corresponds to an F function of the 
type given in Eq. (12) with A^ = \, B = Q and C = 1. 
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Solving the model is straightforwardly carried out using 
the procedure outlined in this article. For any set of 
rapidities which is a solution of the Bethe equations, the 
eigenenergy of the state is then given by 



N 

E = ^A,. (32) 

At zero-coupling, the M-pairs eigenstates are simply 
obtained occupying any M energy levels with a Cooper 
pair (pseudospin up) while the remaining N — M are 
empty. This can be represented by the set of rapidi- 
ties Aq, = which according to Eq. (30) lead to an 
up pseudospin at every energy Cj^ present in the set 
of rapidities. This state obviously leads to diverging 

A(ei^ ) = z2 ^ occupied levels, but lin- 

earizing the Bethe Eqs. (31) tells us that every rapidity 
can be written as = £«. — g at weak coupling. Wc 
consequently solve Eqs. (16) using = gA{ej) which 
are given at zero coupling by 




r 1 occupied ej 
■' y empty ej. ^ ' 

While in other applications, such as nuclear 
physics'^"'''^, a different distribution of the levels could 
be better suited, we present results for equally spaced 
single particle energy levels within a bandwidth D, i.e. 
tj = + ^ This Richardson model is the one typ- 
ically used in condensed matter systems^. Both Aj = 
gA{€j) and the corresponding rapidities for the ground 
state and two excited states are plotted in Fig. 1. 

The black circles show the points at which we actually 
solve the equations, the step size in g between two points 
being Ag/d — ^. These results were computed using the 
6 first derivatives for the regression (see section II). For 
reference we present in gray the numerical solution as a 
continuous function of g/d. The figure therefore shows 
how this approach allows steps in the coupling constant 
which are very large compared to the actual scale on 
which rapidities themselves vary. 

The middle panel shows a state specifically chosen to 
generate a region of rapid variation of the rapidities that 
one can see around g/d = 1.3. When solving directly 
for the rapidities themselves using Eqs. (31), one would 
need an extremely small step size in this region in order to 
maintain stability. In the current approach however, this 
complex structure causes no problem since the variables 
Aj remain smoothly varying functions. This constitutes 
a prime example of the remarkable capacities of this new 
approach which allows the scan in g to be performed 
orders of magnitude faster than using the standard Bethe 
Eqs. (31). 



FIG. 1: gA{ej) (left) and the real part of the corresponding 
rapidities (right) for the ground state (top) and two excited 
states of the equally spaced Richardson model (N=50, M=25). 
The circles are the only points needed when using the 6 first 
derivatives for the regression. 

B. Generalized Dicke model 

We now turn to a second related model, the General- 
ized Dicke model. It describes a collection of N multi- 
level systems coupled uniformly to a single bosonic mode 
and has been shown to be Bethe Ansatz solvablc^'^'^^. 
Here, we again consider the case of 2-level systems (repre- 
sented by a spin 1/2) which is evidently the relevant case 
for cavity-based quantum computing proposals^^. The 
Hamiltonian takes the form 

N N 

where {Sj , S'^ , S~} are the usual generators of spin 1/2 
representations of SU(2) at each site j. The boson fre- 
quency is uj, Ej sets the splitting between both levels of 
every subsystem, while V controls the strength of the in- 
teraction. This Hamiltonian conserves the total number 
of excitations b^b -\- (S'| -|- |). The construction of 
M-excitations eigenstates is achieved using the Gaudin 
creation operators 

S+(A„) = 6t+^^^5+ (35) 

3 ■' 

acting repeatedly on the pseudo- vacuum state |0) = 
|0; I, I, 4-) which contains no boson and is fully down- 
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polarized. The appropriate sets of rapidities {A^}^! 
must fulfill the Bethe equations'^^ 



2^2 2V^ 



N 

E 



1/2 



M 



Here, wc can identify the general form (Eq. 12) with 
g = V^, Ai = 1/2, B = -1 and C = w. The energy 
eigenvalues corresponding to each set of rapidities are 



(37) 



Now we turn to the solution of (36) using the method 
proposed in the previous section. The equivalent equa- 
tions for Aj = V^A{ej) read 

By linearizing the Bethe equations (36), one can show 
that Aq, = Ej — give correct solutions in the limit 

-> 0. These values lead to = w - Cj at ^ = 0. 
As in the Richardson model, Aq, = Cj^ leads to an up- 
spin at level e^^. Another possible solution is obtained 
for \a = ijJ which leads to hj = 0. In this case, every 
\a ^ u) corresponds to an additional bosonic excitation 
at V"^ = 0, since eigcnstates are constructed using the 
operator (35). Any of the d = '^f^i {^) combinations of 
M excited spins and bosons gives us a possible eigenstate 
of i? at = 0, which we deform numerically to find 
eigenstates and eigenvalues at nonzero coupling > 0. 

Figure 2 shows the rapidities and A^- for two different 
states as a function of V^. We made a specific choice 
for the inhomogeneity by using equally spaced splittings 

= — -y + ^ J. This would ultimately correspond to a 
coarse-graining of a flat distribution of splittings within 
a bandwidth D. Moreover, we choose the bosonic fre- 
quency to be at the midpoint of the band, i.e. u = 0. 
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FIG. 2: V'^A{ej) (left) and the real part of the corresponding 
rapidities (right) for the "equally spaced" Generalized Dicke 
model (N=60, M=20). The circles are the only points needed 
when using the 5 first derivatives for the regression. 



H = HoSq + AjSo ■ Sj 



(40) 



j¥0 



While the external magnetic field is only coupled to the 
central spin i = 0, the conservation of the total excitation 
number Sq + ^ Sj makes the model equivalent to 

37^0 



H = hgeS^ + hgNj2s^ 

37^0 



/ ,-'-3 
37^0 



AjSq ■ Sj, 



(41) 



which would also include coupling of the bath spins to 
the external magnetic field. 

In principle, any linear combination of the integrals 
of motion also leads to an integrablc model which com- 
mutes with Hamiltonian (40) and therefore shares the 
same eigenbasis. Consequently any Hamiltonian 



C. Central Spin model 

Finally, we apply the method to the Central Spin (CS) 
model which in its most frequently encountered form de- 
scribes a single spin coupled to an external magnetic 
field and via hypcrfine interaction, to a spin bath. It 
is straightforwardly obtained from the commuting con- 
served quantities of the SU(2) XXX Gaudin models given 
by 

Ri = S^ + Y^^^Si-Sj, (39) 

3^i 

by defining the exactly solvable model as H = with 
the choice eo = 0, /lo = ^ and A, = — j- , i.e. 

9 



(42) 

is also Bethe Ansatz solvable. While this allows the treat- 
ment of models which would include hyperfine interac- 
tions between the bath spins i > 0, integrability restricts 
the allowed values of the couplings. For example, if we 
suppose that bath spins are homogeneously coupled to 
the external field, i.e. rjj = r] Vj > 0, integrability im- 
poses the absence of coupling between bath spins. On 
the other hand, models with non-zero interaction in the 
bath can be built from any choice of rjj but this in turn 
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would force the presence of an inhomogeneous magnetic 
field acting differently on every bath spin. 

In the limit h = 0, one could study models with 
bath interactions provided the couplings respect the con- 
straints 

-^(^,(r?o - Vk) - Akim - Vj)) = AjAk, (43) 

Vj - Vk 

which has important consequences. For example when 
any two nuclear spins are equally coupled (Aj = Ak = A) 
to the central one, the coupling between those two spins 
would also need to be given by Bj^ = A. For a given set 
of distinct Aj , it would however still be possible and inter- 
esting to study Central Spin systems with imiform bath 
couplings Bij = B 01 any integrable case that emerges 
from a given choice of distinct rjj. 

Wc still choose to focus on Hamiltonians of the form 
(40). Is is straightforward to show that it commutes with 
the Richardson Hamiltonian (Eq. 29) defined by param- 
eters g = — ^,eo = 0, ej = — They therefore share 
common eigenstates so that, in principle, solving Eqs. 
(31) for negative values of g would give us the eigenstates 
of the central spin Hamiltonian. However, we follow the 
alternative road of inverting the spin quantization axis 
{z —z) and solving Eqs. (31) for positive g. Eigen- 
states are then be obtained using the operators 

^^w = Ex^' (44) 

acting on the fully up-polarized pseudo- vacuum |0) = 
Itt ■ • ■ t) and the eigenenergy of a given eigenstate is 
given by 

Nt, ^ ^1 h 
iSo({Afc}fc=i,...,M) = Yl T + 2 E \ T + 9 • 

(45) 

In some of the physically relevant systems that can 

be described with Central Spin model (such as Quan- 
tum Dots or Nitrogen- Vacancies (NV) in diamond'^'"'), 
the proper interactions between the CS and the nuclei 
do not give rise to equally spaced values of — 

fact, some of the parameters tj can be very close to one 
another while others are strongly separated. While the 
Bethe equations are exactly the same as for the Richard- 
son model this distribution of has a direct consequence 
on the application of the current method. Since for 
equally spaced levels with separation d, the Bethe equa- 
tions can be written only in terms of gjd, the size of 
the steps in g that one can use is simply controlled by 
the d parameter. However, when a number of levels are 
very close by while others are far from one another, these 
variable spacings lead to variations of A(ej) controlled by 
different scales. Naturally, using a step in g that is a given 



fraction of the smallest distance — would insure 
stability but it would ultimately necessitate a large num- 
ber of points to reach the strong g = \ limit (weak mag- 
netic field). It is therefore beneficial to use a variable step 
size. Starting from small steps at low </, we turn to larger 
ones when g has increased sufficiently and the behavior 
of the solutions vary on a much slower scale. Variables 
steps could also be used in the previously studied models, 
but in the current case it is more or less necessary to do 
so in order to achieve fast computation. 

Figure (3) shows the behavior of the rapidities for an 
NV-Center (iV = 50 and M = 25), in the interval g going 
from to 0.1 the Newton steps were taken for gS{€j) 
with a step size Ag = 1/50, for g > 0.1 the step size 
is increased to Ag = 1/15 {Ag = 1/20 for the excited 
state). The first 6 derivatives are used. 

Due to the structure of the levels Sj we choose to re- 
strict the y-axis in the plots of the rapidities to make the 
intricate structure of the solutions visible. The insets 
present them in the complete range. 




FIG. 3: 3A(ej) (left) and the rapidities (right) for the 
dipole couplings of an NY-Center*® ej = —^/Aj and Aj = 
5.6(^)^(1 - 3cos(6'i)^)[MHz], where a^j = 1.54[A] is the 
nearest neighbor distance for diamond and Rj is the distance 
between the j-th carbon atom and the defect (center) . 6j is 
the angle between the magnetic field and Rj. The vectors 
Rj are corrected by a small amount of randomness to avoid 
degeneracies. Top panel shows the ground state; bottom one 
shows an excited state. Insets show the full range of rapidi- 
ties while the main figures focus on the region with non-trivial 
structure in the given range of magnetic field. 

Once again, the black dots represent the points at 
which the equations are solved. One can notice the in- 
crease in the step size at large g = j^. In spite of the 
necessity of using smaller steps, these plots show with- 
out a doubt that this approach still allows \is to realize 
the scan in g using large steps on the scale on which 
rapidities themselves vary. 

A variable step-size could naturally be defined by mak- 
ing use of the fact that we compute the truncated Taylor 
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expansion up to a given order. One could fix the appro- 
priate steps simply by defining a desired level of precision 
for the limited expansion. This would allow one to con- 
trol stability, while adapting the steps to the behavior of 
the A{€j) variables at the current g point. 

V. CONCLUSIONS 

In this work we have shown how, for rational Gaudin 
models, the set of non-linear coupled Bethe equations 
can be transformed to a new set of equations which is 
significantly simpler to solve numerically. We presented 
the complete algorithm and applied it to a variety of 
models derived from the generalized Gaudin algebra. 

This work exclusively presents and demonstrates the 
capacities of this new approach to the solving of Bethe 
equations. In light of previous work^^"^^, it should how- 
ever be clear that by allowing rapid and systematic solv- 
ing of a decent fraction of the full Hilbert this can directly 
be used for relevant physical calculations. For example. 



studies of the non-equilibrium dynamics of the discussed 
models are currently pursued by the authors. The out- 
come of these calculations should prove invaluable con- 
sidering that this technique gives direct access to regimes 
which have been hard to describe before. To name only 
this one, the weak magnetic field limit of the Central spin 
models is a prime example. 

While we focused on rational Gaudin models, algo- 
rithms tailored to the trigonometric or hyperbolic XXZ^" 
models or general XYZ^'' could possibly be built along 
similar lines. This is a more fundamental question that 
is left open for the time being. 
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